Coherent phonon manipulation in coupled mechanical resonators 
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Coupled mechanical oscillations were first 
observed in paired pendulum clocks in the 
mid-seventeenth century and were extensively 
studied for their novel sympathetic oscillation 
dynamics[l|, Q|. In this era of nanotechnologies, 
coupled oscillations have again emerged as sub- 
jects of interest when realized in nanomechani- 
cal resonators for both practical applications and 
fundamental studies [3 11]. However, a key ob- 
stacle to the further development of this archi- 
tecture is the ability to coherently manipulate 
the coupled oscillations. This limitation arises as 
a consequence of the usually weak coupling be- 
tween the constituent nanomechanical elements. 
Here, we report parametrically coupled mechani- 
cal resonators in which the coupling strength can 
be dynamically adjusted by modulating (pump- 
ing) the stress in the mechanical elements via 
a piezoelectric transducer. The parametric con- 
trol enables the coupling rate between the two 
resonators to be made so strong that it ex- 
ceeds their intrinsic energy dissipation rate by 
more than a factor of four. This ultra-strong 
coupling can be exploited to coherently trans- 
fer phono n p opulations, namely phonon Rabi 
oscillations[12j, 1 1 31 ] . between the mechanical res- 
onators via two coupled vibration modes, real- 
izing superposition states of the two modes and 
their time-domain control. More unexpectedly, 
the nature of the parametric coupling can also 
be tuned from a linear first-order interaction 
to a non-linear higher-order process in which 
more than one pump phonon mediates the co- 
herent oscillations. This demonstration of multi- 
pump phonon mixing echoes multi-wave photon 
mixing [14] and suggests that concepts from non- 
linear optics can also be applied to mechani- 
cal systems. Ultimately, the parametric pump- 
ing is not only useful for controlling classical 
oscillations[15] but can also be extended to the 
quantum regime[12, 13, ltH- 18|. opening up the 
prospect of entangling two distinct macroscopic 
mechanical ob jects [191 l20l|. 

The dynamic parametric coupling is developed in 
GaAs-based paired mechanical beams shown in Fig. QJi, 
in which the piezoelectric effect is exploited to medi- 
ate all-electrical displacement transduction 2l( . The fre- 



quency response of beam 1 measured by harmonically 
driving it while the parametric pump is deactivated dis- 
plays two coupled vibration modes (Fig. Q~p) , where mode 
1 (o>i = 2tt x 293.93 kHz) is dominated by the vibra- 
tion of beam 1 while mode 2 (ui2 = 2ir x 294.37 kHz) 
is dominated by the vibration of beam 2. The ampli- 
tude of mode 2 is much smaller than that of mode 1 
reflecting the energy exchange due to the structural cou- 
pling via the overhang is inefficient because of the eigen- 
frequency difference between the two beams. This differ- 
ence can be compensated by activating the parametric 
pump, which is induced by piezoelectrically modulating 
the spring constant of beam 1 with the pump frequency 
Lu p at around the frequency difference between the two 
modes, Alu = uj 2 — ^i (Fig. QJ). 

The dynamics of this system can then be expressed by 
the following equations of motion: 

Xi + 7iii + [w 2 + Ti cos(u) p t)]xi 

+A cos{uj p t)x2 = Fx cos(uidt + (j)) (la) 

'±2 + 72^2 + [W 2 + T 2 COs(0J p t)]x2 

+Acos(uj p t)xi = F 2 cos(uidt + (j)), (lb) 

where Xi [i = 1, 2) is the displacement of the i-th mode, 
uji is the mode frequency, ji is the energy dissipation 
rate, Fi is the drive force (iq ^> F^), and iOd is the drive 
frequency. When the frequency mismatch between mode 

1 and mode 2 is compensated by activating the pump at 
ujp ~ Aw, the terms containing A transfer phonons (os- 
cillations) from one mode to the other (Fig. [Hi). This 
inter-modal coupling can also be explained by the mixing 
of mode 1 (2) and the Stokes sideband, 0J2 — U! p (the anti- 
Stokes sideband, u>i + uj p ) of mode 2 (1) leading to the 
normal-mode splitting in the strong-coupling regime (2lT- 
I25I ] (Fig. [TJi) . The equations also include terms propor- 
tional to r^. These terms lead to intra-modal coupling 
which becomes significant for the higher-order paramet- 
ric coupling as shown later. The above model can re- 
produce all the experimental results and is described in 
detail in Supplementary Information. 

The parametric normal-mode splitting for modes 1 and 

2 was experimentally confirmed by applying a pump with 
voltage V p at u> p in addition to the weak harmonic driving 
at uid to beam 1 (Fig. [1^) . Probing the modes via beam 1 
while the parametric pump is activated shows that mode 
1 splits into two when ui p ~ Au = 2tt x 0.44 kHz clearly 
demonstrating the strong dynamical coupling (Figs. [2jt- 
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FIG. 1: Paired mechanical resonators and the parametric pumping scheme, a, Schematic drawing of the sample, the 
measurement and the false-colour scanning electron micrograph. The two doubly clamped beams are structurally interconnected 
via the coupling overhang, b, Frequency response of beam 1 measured by driving it with Vd = 1.5 mV p _ p while the pump 
is deactivated (V p = V p _ p ). The mode shapes at uj\ and LJ2 obtained by finite element method calculations are also shown 
(exaggerated for clarity), c, Schematic drawing of the parametric pumping protocol in an equivalent spring model, ki, fe, 
and k c are respectively the spring constant of beam 1, beam 2, and the coupling overhang. F represents the parametric pump 
amplitude, d, Schematic drawing of the parametric pumping protocol in an energy diagram. 



[5p). The phonon transfer from mode 1 to mode 2 can be 
confirmed in the response of beam 2 that was detected at 
frequency Wd + w p (Figs. [2jl-[2f). The results indicate the 
emergence of its vibration and mode splitting at w,j ~ wi 
and uip ~ Aw showing that mode 2 (u>2 = Wi + Aw) is 
excited by the dynamical coupling. A phonon reaction 
picture can help to understand this elementary process 
where phonons are created in mode 2 at the expense of 
probe phonons in mode 1 and pump phonons, i.e., via the 
one-pump phonon absorption process, Hwi + ftio p — > fiw2 
(Figs. [TJi and[3h). Together with the reversed emission 
process, TibJi — > tkoi + huj pi strong coupling is established 

(Fig. mo. 

The V p dependence of the mode splitting at w p = Aw 
shows that the coupling strength is highly controllable 
with the splitting being proportional to V p (Fig. HJ). 
This linear V p dependence is due to the fact that the 
inter-modal coupling coefficient, A, is proportional to V p , 
which can be theoretically reproduced by Eqs. [Ta] and 
[Tbl (Fig. [2j<) . The separation between the split peaks 
provides the coupling rate, g, which can become strong 
enough to exceed the intrinsic energy dissipation rate of 
the two modes (71 ~ 72 = 2tt x 22 Hz) by more than a 
factor of four (g — 2ir x 90 Hz for V p = 1.0 V p - P ), i.e., 
ultra-strong parametric coupling can be achieved in this 



system. 

More remarkably additional mode splitting, in which 
the pump frequency does not correspond to the frequency 
difference between the two modes, can also be observed. 
The additional splitting occurs when w p ~ Aw/2 = 
2tt x 0.22 kHz for both modes 1 and 2 (Fig. Gfc). This 
splitting is caused by a second-order parametric coupling 
via two-pump phonon absorption/emission process, i.e., 
fiwi + 2fiw p thu)2, which leads to dynamic coupling 
between mode 1 (2) and the second Stokes sideband, 
w 2 — 2w p (the second anti-Stokes sideband, Wi + 2w p ) of 
mode 2 (1). The V p dependence of this mode splitting in- 
dicates that it has a parabolic dependence (Fig. (2). This 
is because the second-order process requires a two-step 
phonon excitation path from mode 1 to mode 2, and vice 
versa, through the intermediary energy level (wi +W2)/2 
via both the intra- modal coupling (r, oc V p ) and the 
inter- modal coupling (A oc V p ), i.e., [Ti x A] oc VI 2 , as 
shown in Fig. [3j. Again, the corresponding mode split- 
ting shows good agreement with the theoretical simu- 
lations (Fig. 0). Uniquely, the second-order phonon 
processes can be observed in the present system owing 
to the strong parametric coupling between the two res- 
onators/modes. 

The ultra-strong parametric coupling in the paired me- 
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FIG. 2: Normal-mode splitting induced by parametric pumping, a-c, The drive frequency (wd) and the pump frequency 
(cj p ) response of beam 1 detected at frequency uJd for three different pump voltages, V p — 0, 0.5 and 1.0 V p _ p . d-f, The u>d 
and aip response of beam 2 detected at frequency uid + u p for V p = 0, 0.5 and 1.0 V p _ p . g,h, Simulation results for (c) and 
(f), which were performed for the theoretical model expressed by Eqs. Hal and [TBI i,j, The V p dependence of the splitting of 
mode 1 induced by the first- and second-order parametric coupling at uj p = Auj and Aw/2, respectively. The broken curves 
in (i) and (j) represent the theoretical values of the mode splitting, which are proportional to V p and (Vp) 2 respectively (see 
Supplementary Information). The mode splitting corresponds to the coupling rate, g. k,l, Simulation results for the splitting 
of mode 1 for uj p — Au and Aoj/2, respectively. 



chanical resonators opens up a path to coherent control 
of the mechanical oscillations. The time-domain mea- 
surements utilizing the pulse sequence shown in Fig. [3^, 
and described in Methods Summary enable us to observe 
coherent and periodic energy exchange between the two 
beams/modes. The pump frequency dependence of the 
time-domain response of beam 2 at u>2 clearly shows the 
periodic amplitude oscillations at lo p ~ Alj for the first- 
order (n = 1) parametric coupling (Fig. [3p). The V p 
dependence at uj p = Aw shows that the vibration energy 



of mode 1 (beam 1) can be transferred to mode 2 (beam 
2) and back eight times at V p = 1.0 V p _ p (Fig. |3Ji). 
Coherent energy exchange for the second-order (n = 2) 
parametric coupling can also be observed at oj p ~ Aw/2 
(Fig. [3p), where up to five oscillation periods are ob- 
served in the range of V p < 1.0 V p _ p (Fig. [3b). The 
coupling rate g, extracted from the Fourier transforms of 
the time-domain response is proportional to V p for the 
first-order coupling and (V p ) 2 for the second-order cou- 
pling (Fig. [4J. These coupling rates correspond perfectly 
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FIG. 3: Coherent energy exchange between the two beams/modes, a, Pulse sequence used for the measurements. b,c, 
The u>p dependence of the time-domain response of beam 2 at CJ2 when the system is subjected to the above pulse sequence, 
from the experiment (b) and the corresponding simulation (c) for V p = 0.5 V p _ p . d-g, The V p dependence of the time-domain 
response of beam 2 at 0J2 for u v — Acj/n, where n = 1, 2, 3 and 4, respectively. The higher-frequency fringe patterns found 
for n > 2 are caused by the overlap of the coherent oscillations for the (n-l)-th order coupling, which is confirmed in (b). 
h-k, Schematic drawings of one, two, three, and four-pump phonon absorption processes, where A is the inter-modal coupling 
coefficient while Ti and F2 are the infra-modal coupling coefficient for mode 1 and mode 2, respectively. 



to the mode splitting observed in the frequency response 
measurement. 

Even more surprisingly, the time-domain measure- 
ments enable us to observe higher-order parametric cou- 
pling, i.e., n > 3. Figures [3f and [3g show the coher- 
ent energy exchange between the two beams/modes for 
lj p = Acj/3 = 2ttx 0.147 kHz and uj p = Aw/4 = 2ttx0.11 
kHz, respectively. These coherent oscillations are caused 
by the n-pump phonon absorption/emission processes, 
i.e. huji + nHujp ■<-> Hu)2, through the intermediary en- 



ergy levels via intra- and inter-modal coupling, e.g., 
[ri x Ti x A], as shown in Figs. [3] and[3]t. The Fourier 
transforms of the time-domain response reveal that the 
coupling rate exhibits (V p ) n dependence even for n > 3 
(Fig. fJJ, which again shows good agreement with the 
theoretical model. 

The present results show that electromagnetic pulse 
techniques, which are commonly used to coherently ma- 
nipulate quantum two- level systems [l^. [2fj| . can also 
be applied to coherently control mechanical systems, 
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FIG. 4: Coupling rate for the n-th order parametric 
pumping. The logarithmic plots of the coupling rate, g, for 
u) p — Auj/n (n = 1, 2, 3 and 4) with respect to V p , ob- 
tained from the Fourier transforms of the time-domain re- 
sponse shown in Figs. [3ji-[3g. Solid lines represent the the- 
oretical values, which are proportional to (V p ) n (see Supple- 
mentary Information). 



which enables the realization of superposition states be- 
tween the two coupled vibration modes. By tuning the 
parametric-pump frequency, multi-wave phonon mixing 
involving an arbitrary number of pump phonons can be 
achieved in an analogous fashion to multi-wave photon 
mixing 14]. The parametric pumping allows highly con- 
trollable time-domain manipulation of the phonon pop- 
ulation with the adjustment of the pump-pulse dura- 
tion, thus permitting tt and 7r/2-pulse operations on the 
Bloch sphere [IH (see also independent experiments at 
LMU with a single mechanical resonator 27]). This co- 
herent control further expands the applications of me- 
chanical resonators, including the high-speed operation 
of high-Q mechanical resonators [15j and mechanical logic 
operations [3]. Although the system demonstrated here 
is in the classical regime with large mode occupation, 
these techniques can also be extended to the quantum 



regime 12 , 13, HI > leading to the exciting possibility of 
quantum-coherent coupling and entanglement between 
two distinct macroscopic mechanical ob i ects flil. I20L [2rj| . 



METHODS SUMMARY 

Fabrication. The sample was fabricated by pho- 
tolithography from a heterostructure consisting of 
Alo.25Gao.75As, Si-doped n-GaAs, undoped i-GaAs and 
2-/Ltm-thick Alo.65Gao.35As sacrificial layers grown on a 
GaAs(OOl) substrate by molecular beam epitaxy (Fig. 
Hk). AuGeNi was deposited on the supporting part to 
obtain an ohmic contact to the conductive n-GaAs layer, 
while Au/Cr gates were formed on the top of the beams 



(Fig. QJi). The suspended structure was completed by 
deep mesa and isotropic sacrificial layer etching, where 
the 40-/im-separated beams were electrically isolated by 
the shallow mesa etch. All the measurements were done 
by setting the sample in a vacuum (5 x 1CP 5 Pa) at 1.5 
K with a 4 He cryostat. 

Frequency-domain measurement. Beam 1 was 
driven by a continuous AC voltage (Vd) applied to gate 
D with the frequency tod ~ <jJ\- Parametric pumping 
was achieved by applying a continuous AC pump volt- 
age (V p ) to gate P with the frequency lo p ~ 0J2 — ui\. 
The responses of beam 1 and beam 2 were simultane- 
ously monitored via gates Ml and M2, where the signal 
was amplified with a room-temperature low-noise pre- 
amplifier (NF:SA-220F5) and detected with a lock-in am- 
plifier (SRS:SR844) (see Supplementary Information). 

Time-domain measurement. A sinusoidal drive 
pulse voltage (Vd) with cud = uii was applied to gate 
D with the pulse period of 0.1 s. A sinusoidal pump 
pulse with uj p = (cl>2 — wi)/n was applied to gate P at 
the moment the drive pulse was switched off, where the 
pump pulse also has the pulse period of 0.1 s. The time 
response of beam 2 at ll>2 was measured through gate 
M2 with an oscilloscope (Agilent:DSO90254A) via the 
pre-amplifier and the lock-in detector (see Supplemen- 
tary Information) . The data was averaged 20 times with 
the oscilloscope. 
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Supplementary Information for 
Coherent phonon manipulation in coupled mechanical resonators 



MEASUREMENT SCHEME 



Piczoelectrically frequency tunable paired mechanical resonators were prepared to demonstrate parametric coupling 
(Fig. IS 1[) . Beam 1 and beam 2 are structurally (elastically) coupled via overhang parts (see Fig. la in the main 
text) that were formed by the sacrificial layer etching process. The structural-coupling constant can be estimated 
from the results of the DC detuning experiments (described in the next section), where gate T is used to change the 
eigenfrequency of beam 2 with the application of DC voltage (Vf) using a function generator (NF:WF1974) (Fig. IS 1 1) . 
For the parametric pumping measurements, both the drive voltage and pump voltage were applied to beam 1 with the 
function generator. A weak AC drive voltage (Vg — 1.5 mV p _ p ) was applied to gate Dl while the conductive n-GaAs 
layer (see Fig. la in the main text) was grounded via the ohmic contact formed on the supporting part (not shown). 
An AC pump voltage (V^,) was applied to gate P, where the pump frequency (lu p ) was set at around the frequency 
difference between the two coupled vibration modes (u)2 — uii), while Vt was set to V. In the frequency-domain 
measurements, the responses of beam 1 and beam 2 were simultaneously monitored via strain-induced piezoelectric 
voltage at gate Ml and gate M2 with low-noise pre-amplifiers (NF:SA-220F5) and lock-in amplifiers (SRS:SR844). 
Here, we used two lock-in amplifiers to measure each beam's response, one to measure the signal at the drive frequency 
ojd and the other to measure the idler at the frequency u>d + u p as shown in Fig. IS 11 In contrast, in the time-domain 
measurements, the real-time response of beam 2 at the idler frequency Ud + uj p was monitored with an oscilloscope 
(Agilent:DSO90254A) via a lock-in amplifier as shown in Fig. IS11 



STRUCTURAL COUPLING 



In the paired mechanical resonators, the degree of vibration coupling strongly depends on the eigenfrequency 
difference between the two resonators. This was confirmed by applying DC detuning voltage (Vt) to beam 2 while 
beam 1 was harmonically driven (Fig. IS1|) . The Vt dependence of the frequency response of beam 1 and beam 2 for 
V p = Vp_ p is shown in Figs. IS2b andlS2b. respectively. At Vt — V, the amplitude of mode 1 (oji — 2ttx 293.93 kHz) 
is much larger than that of mode 2 (w 2 = 2ir x 294.37 kHz) in beam 1 (Fig. IS2a h This is because the eigenfrequency 
of beam 2 (0 2 ) does not coincide with that of beam 1 (Oi) at Vt = V. The separation of the two mode frequencies 
reduces with the application of the negative Vt and it is minimized at Vt — —0.5 V, showing the avoided crossing of 
the two modes (Figs. IS2h and IS2b). This is because O2, which is initially higher than Oj owing to the fabrication 
error, decreases with increasing negative Vt due to the piezoelectric effect and perfect tuning (f2 2 = Oi) is achieved at 
Vt = —0.5 V. In the perfectly tuned condition, the amplitudes of the two modes (and the two beams) coincide with 
each other, where the lower- frequency mode (w s = 27r x 293.85 kHz) corresponds to the symmetric vibration and the 
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FIG. SI: A topographic optical image of the paired mechanical resonators and a schematic of the circuit for the measurements. 
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FIG. S2: DC detuning voltage (Vt) dependence of the frequency response of (a) beam 1 and (b) beam 2 measured by harmonically 
driving beam 1 with no parametric pumping. Dotted lines show the eigenfrequency of beam 1 (Sli ) and beam 2 (fb). The 
response at Vt = V is highlighted by a broken line. Drawings of the symmetric vibration and anti-symmetric vibration for 
the perfect tuning point are also shown. 



higher-frequency mode (u a = 2ir x 294.08 kHz) corresponds to the anti-symmetric vibration. The difference between 
the two mode frequencies at the perfect tuning point provides information about the structural coupling constant 
(c = n 2 c ), where Q, c = ^(uj 2 a -uj 2 )/2 = 2ir x 8.22 kHz is 2.8 % of From Fig. [52ji (and EH)) the piezoelectric 
frequency detuning coefficient is estimated to be 5Vl 2 /8Vt = 2ttx 0.75 kHz/V. 



THEORETICAL SIMULATIONS 



Model 



We used a standard model of two coupled harmonic oscillators with the linearly voltage-controlled detuning of 
resonance frequency. The equations of motion are given by 

X1+71X1 + p +aV 1 (t)]X 1 = c{X 2 - + Fi(t) (Sla) 
X 2 + l2 X 2 + p 2 + aV 2 (t)}X 2 = cpfj - X 2 ) + F 2 {t). (Sib) 

Here, Xi (i — 1,2) is the displacement of the i-th oscillator, Hi and ft 2 = O x + Sfl (Qi >> 8SY) are the resonance 
frequencies of beam 1 and beam 2, respectively, 7^ = Qi/Qi is the energy dissipation (damping) rate where Qi is the 
quality factor, V$ is the detuning gate voltage, Fi is the harmonic driving force, and c is the geometrical coupling 
constant between the two oscillators. The term aVi(t) gives the eigenfrequency detuning and the device we used has 
the detuning coefficient of a/2fl; = Sfti/SV = 2irx 0.75 kHz/V (Figs. \S2h andlS2b). We describe the case when beam 
1 is harmonically driven at frequency cud with amplitude Fi — F, the eigenfrequency of beam 1 (fii) is parametrically 
modulated (pumped) at frequency uj p with a voltage amplitude of V p — T/a, and no detuning gate voltage is applied 
to beam 2. Then the equations of motion become 

Xi +7X1 + pl+Tcos(u p t)]Xi = c(X 2 - Xi) + F cos(a; rf t + <f>) (S2a) 
X 2 + 7X2 + [ill + 2n 1 Sil}X 2 = c(Xi - X 2 ), (S2b) 

where <fi is the phase factor determined from the phase relation between the harmonic driving F cos^^t) and the 
parametric pump Tcos(u! p t), and the damping rate was assumed to be 7 = 71 =72- When no parametric pump is 
applied, the equations become 

X l + 1 X l + (p + c)Xx -cX 2 = F cos(uj d t + <P) (S3a) 
X 2 + 1 X 2 + (p + 2fti«] + c)X 2 - cXi = 0, (S3b) 
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d^ + ^dt + ^ 



x 1 
x 2 



F cos(uJdt 




x 2 



Using the diagonalization orthogonal matrix U dchncd by 

1 /VXTS x/A - d 



U 



(S4) 



(S5) 



where d = £liS£l and A = ^ c 2 + n 2 5£l 2 . Eq. [S5]is diagonalized as 



rf2 d o2 



■f'i 

2:2 



A 

-A 



Xi 



U 2 i 



F cos(ojdt 



where Xi is the displacement of the z-th mode given by 



X\ 
J'2 



= u 



x 2 



The equivalent expression of Eq. IS6I is 



x\ + 7±i + uj 1 x 1 



x 2 + 1x2 + w 2 x 2 



Ui\F cos(uidt 
U21F cos(ui d t 



= Ql + c + d- Vc 2 + d 2 
= o 2 + c + d + Vc 2 + d 2 . 



(S6) 



(S7) 



(S8a) 
(S8b) 



The two independent modes are referred as mode 1 and mode 2 hereafter. 

In the case when the parametric pumping is ON, Eqs. IS2al and IS2bl can be expressed by using the two modes as 



xi + 7x1 + uifxi + (T1X1 + Ax 2 ) cos(u! p t) = Ux\F cos(uidt + <j>) 
X2 + 7±2 + u\x 2 + {^2X2 + Aati) cos(u} p t) = U 2 iF cos(wdt + 0), 



(S9a) 
(S9b) 



where Y\, T 2 , and A are given by 



Ti = rut, = 



r 2 = rctf, = - 1 



a = ru n u 21 = 



sjc 2 + n 2 sn 2 1 

riiSfi 
.jc 2 + n 2 5n 2 \ 

Yc 



2Vc 2 + njsn 2 ' 

By simplifying Eqs. IS9al and IS9bl with F\ = U\\F and F2 = U21F, we obtain the forms of Eq. 1 in the main text, 



x\ + 7.T1 + \ui\ + Ti cos(w p t)]xi + Acos(w p i)a: 2 = F\ cos(uidt + 1 
X2 + J&2 + [w 2 + r 2 cos(a;pi)]2: 2 + Acos(wpi)a;i = F 2 cos(uid.t + > 



(la) 
(lb) 



The terms including A represent the inter-modal coupling, which transfers phonons from one mode to the other. On the 
other hand, the terms containing I\ represent the intra-modal coupling, which becomes significant for the higher-order 
parametric coupling as described in the latter section. Because Y = aV p , both the intra-modal coupling coefficient Yi 
and the inter-modal coupling coefficient A are proportional to the pump voltage V p . Note that A is maximized when 
50 = but still sufficiently large when <50 is comparable to c/fii, where c ~ 2.6 x 10 9 Hz 2 , 51i ~ 1.8 x 10 6 Hz, and 
SD, = 2.5 x 10 3 Hz in the present paired mechanical resonators; therefore, <5£! ~ c/fii. 
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G 
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FIG. S3: Numerical simulations of oscillation amplitudes of [(a)-(c)] beam-1 at drive frequency u)d, [(d)-(f)] beam-2 at the 
first-order idler frequency u>d + and [(g)-(i)] beam-2 at the second-order idler frequency ujd + 2oj p as a function of pump and 
drive frequency. 



Higher-idler rotating-frame approximation 



Equations la and lb can be numerically solved by using rotating-frame approximation with higher-harmonic idler 
components. By introducing slowly varying m-th harmonic idler amplitudes, a m (t) and b m (t) with m as any integer, 
x\ (t) and X2 it) can be expanded as 



x\{t) = Re 
X2(t) = Re 



+oo 



^2 %(*)e i(w+m "' )i 

.= — OO 

b m (t)e l{u+m ^ )t 



(Slla) 
(Sllb) 
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FIG. S4: Simulated pump frequency dependence of mode 1 by (a) excluding the inter-modal coupling (I\ 7^ 0, A = 0), (b) 
excluding the intra-modal coupling (r, = 0, A 7^ 0), and (c) including both intra- and inter-modal coupling (r, 7^ 0, A 7^ 0). 



Then, the equations of motion become a series of equations: 



[2i(u d + muj p ) + 7] a m + [-(u d + muj p ) 2 + i-i{u d + mu) p ) + uf\ a Tl 



[2i(u> d 



Y\ ]?i A A 

l]b m + [-{uJd + muip) 2 + i~f{u)d + muj p ) + w|] b m 

r 2 r 2 a a 

+ -^b m+ i + — o m _x + — a m +i + -^ a m-i — 0, 



(S12a) 



(S12b) 



where the second-order time derivatives of a m and b m were neglected. The solution in iV-th order approximation can 
be obtained by solving the equations by neglecting the high-order idler, a m and b m for m > N. 



Simulations of frequency response and idler spectrum 



The frequency response and the idler spectrum, i.e., the time- independent solutions, can be obtained by neglecting 
the first terms on the left-hand side of Eqs. IS12al and IS12b[ The solution for AN + 2 vibration amplitudes a m ,b m 
can be obtained by calculating the inverse of the (4iV + 2) x (4AT + 2) matrix. The details of the calculation will 
be reported elsewhere and we here only show the results of ninth-order approximation (N — 9) for the frequency 
response of beam 1 (U^ao + U^fbo) and the first- and second-order idler amplitude of beam 2 (U^a+i + U^b+i 
and t/ 2 ^ 1 a+2 + with three different pump intensities in Figs. IS3a - IS3t . Here, we introduced the unit pump 

intensity To = 2ttuj 2 /Q. The condition of strong coupling is given by T > Tq. The avoided crossing was observed 
in the frequency response when u> p = 102 — oj\ (Fig. IS3b ). This coupling results in phonons being created in mode 
2 [102 ) at the expense of probe phonons in mode 1 (wi) and the pump phonons (w p ), i.e., via the one-pump phonon 
absorption process, Huix + fvjj p —¥ Tujji- The creation of phonons at o>2 is confirmed by the first-order idler amplitude, 
which appears when u) d — Wi and uj p = u>2 — u)\ (Figs. IS3t s and lS3f ). On the other hand, the second-order parametric 
coupling via the two-pump phonon absorption process, tkJ\ + 2frio p — > fvjj2, appears when uj p = {u)<2 — Wi)/2. This 
process is confirmed by the second-order idler amplitude, which appears when Ud = U3\ and uj p = (lu2 — w±)/2 (Figs. 
EShandEl)- 



Inter- and intra-modal coupling 



The role of the inter- and intra-modal coupling can be understood by the simulation for a model, which excludes 
the inter- or intra-modal coupling terms. When the inter-modal coupling terms are neglected, i.e., 7^ 0, A = 0, 
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normal-mode splitting does not occur at uj p ~ 0J2 — (Fig- IS4b) . Instead, ± nw p lines appear in the low pump 
frequency region due to the self-modulation of the i-th mode frequency with u p , i.e., I\ cos(u> p t)xi. In contrast, when 
the intra-modal coupling terms are neglected while the inter-modal coupling is taken into account, i.e., I\ = 0, A 7^ 0, 
normal-mode splitting occurs at co p ~ 0J2 — Wi, whereas uji ± nw p lines disappear (Fig. IS4b ). This reveals that the 
inter-modal coupling terms including A plays an essential role in the parametric coupling. It should be noted that 
the normal- mode splitting appears only at lj p ~ U2 — u)\ and not at to p ~ (u)2 — Wi)/2 in this model excluding the 
intra-modal coupling terms. The normal-mode splitting at uj p ~ (0J2 — w\)/2 appears only when both the inter- and 
intra-modal coupling terms are taken into account (Fig. IS4b ). indicating that the intra-modal coupling is necessary 
for the higher-order coupling. 



Analytical formula for mode splitting 

We derive the analytical formula for mode splitting as a function of pump voltage. We start again from the equations 
of motion, Eqs. IS9al and IS9bl and perform the perturbation expansion with respect to the pump amplitude, Ti, I^, 
and A: 

xi + 7x1 + io\x\ + (T1X1 + A1C2) cos(uj p t) — UuF cos(u>dt + <f>) 
it 2 + 7±2 + <^2 X 2 + (F 2 X2 + Axi) cos(u) p t) — U 2 iF cos(u!dt + (/)). 
The above equations can be rewritten as 

X l\ 771 / J. , ±\ ( Ui 

where 



(Do + D p ) f 2 J = F cosfat + ft I ) , (S13) 



D = (w+^Tt+^i 1 r» = Z' Ti cos(w p t) Acos(w p i) 

° ^ ^+ 7 | +w ay ' " p ^ Acos(^) r 2 cos(^) 

The formal solution for this equation is given by 

Xl ) = (D o + D p )- 1 Fcos(^ + 0) ( ^ n 

x 2 J \ U 2 \ 



m=0 

The operator D _1 is defined as 



(Do -1 - Do^DpDcT 1 + Do^DpDo^DpDo 1 ) Fcos{uj d t 

U21 



£ (-ifDo- 1 (D p D - 1 ) m ^cos(^ + ^) ( ^ n V (S14) 



where Xi( a - 1 ) an d X2 (<•<->) are the mechanical susceptibilities of the two modes: 

1 \- 1 1 r a- 1 1 



(Tl(oj) — CJ 2 + + ' CT2(^) — + 170; + W| 

Using Fourier transformation, 



ii(w) \ = [°° ( x i( t ) 



then 



flu.- 1 = / FvtKU',,1 +<t>) ^ j^j ] .■-'-•'<// 

oi<t>JPX(,,, _ ,,, /> f 1/11 



7re^F(S(w - [ ^ ) + (w d -> -a*, -> -<£) (S16) 
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(a) 



{b } 

I 



U, ](w d J > * 



( d ) -A«V/2 -AHV/2 -A(<a>)/2 

| J, | 



"n(»d) ► • > 



FIG. S5: Schematic diagram of the (a) Oth-, (b) 1st-, (c) 2nd-, and (d) 3rd-order approximation expressed by Eqs. IS2HS241 
respectively. 



2/M \ = ( Xl(u)v{u) 



Do 1 ^ = « " ? " (S17) 



The solution can be expanded in the power series of pump operator D p as 

= D *f (w) - D 1 D P D J f («) + D- 1 D p D- 1 D p D- 1 f(o;) + • • • . (S19) 

The first term provides the 0-th order approximation. For simplicity, in the following, <f> is assumed to be zero. The 
unpumped solution is then given by 

D„ l f(u>) = nFS(u; Ud ) ( ) + K -> ~"d)- (S20) 

The solution has a nonzero value only at = ±w<i, corresponding to the signal oscillation, and the Lorentz-shape 
frequency response of the amplitude is given by the susceptibility, Xi( w d) or X2(^d)- F° r ~ wi and UJ2 — tJi >> 7, 
the contribution from X2(w) is negligibly small and only the lower- frequency mode is excited. We can represent this 
as the diagram shown in Fig. IS5a and approximately obtain 

\-lci. A -TPSit. . , . \„. i, . \ ( Un 



D„ A f(w) = nFS(u> - u d )xi{oJ d ) I ^ j + (u d ^ -u d ). (S21) 
In a similar way, the dominant contribution from the first-order terms is represented as the diagram shown in Fig. 

D " 1 D p D " 1 f(a;) = kFS(uj -u) d - ^ P )x2(^d + u p )^xi(^d) ( ^ + (u>d -> ~^d). (S22) 
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FIG. S6: Schematic drawing of the (a) one-pump phonon absorption and emission processes and (b) two-pump phonon absorp- 
tion and emission processes. 



The other term is at off-resonance and the susceptibility is negligibly small so that the term was neglected. The 
in-coming blue arrow in Fig. IS5b shows that drive frequency Ud is up-converted to u d + u p and has a resonance at 
mode 2. The dominant contribution from the second-order terms is represented as the diagram shown in Fig. 

Do 1 DpDp 1 D p Do 1 f(o;) = nFS(u - u d )xi(^d)^X2(^d + u P )^Xi{ud) ( ^ 

+{u d -> -u d ). (S23) 
Similarly, the dominant contribution from the third-order terms is represented as the diagram shown in Fig. IS5H , 
Do^pD-^pD-^pD- 1 !^) = 

ttFS{uj - uj d - uj. p )x2{ud + u p )^xi(^d)^X2{u d + Up)^xi(ud) (u n ^) 

+ {u d -> -u d ). (S24) 

By summing up all the contributions, we finally obtain 

1 = ttFS(u - U d )- 1 



x 2 {u) J " '"' >l - xi(u d )Ax2(u d + u p )A/4 V 

771 r / \ X2(u d + u p )Axi{u d ) ( \ , , , , 

- TrFS(u-Ud-u p )- — -— - — - — rrj7\ n )+(^d^-u d ). S25) 

We now focus on the first term, which gives the renormalized (dressed) spectrum of mode 1. The second term of 
its denominator is the common ratio of the geometric series and diagrammatically represented as in Fig. IS6k . For 
Up = U2 — U\ and u p » Ud — u\ » 7, we obtain 

Xii^dWn 

(u d -> -u d ) 



nFS(u 


- u d 


ttF8(u 


- u d 


ttFS(u 


- u d 


nF8(u 


- u d 



1 - xiMAx 2 (wd + u p )K/A 

a 2 (ud + Up)Ui 1 ^ 

ai{ud)cr 2 {ud + Up) - A 2 /4 
— 2u2SuUn 



(—2ui5u + i-yui)(— 2u 2 Su + 271^2) — A 2 /4 
— 2u28uU\i 



(u d -t -u d ) 

+ (u d -> -u d ) 



Auiu 2 (5u) 2 — 2ijuiU2Su — A 2 /4 



+ (u d -> -u d ). (S26) 



In the strong-coupling regime (j/5u —¥ 0), it has a maximum when the denominator becomes minimum, i.e., the 
mode splitting 2<5u/ 1 \ which corresponds to the coupling rate, g tyl \ is given by 
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In the case of two-pump phonon process, i.e., uj p — (0J2 — Wi)/2, off-resonance intermediate states at oj = cod + uj p = 
+ ( ^2 — ) /2 should also be taken into account . The diagrammatic representation of the common ratio of geometric 
series is given by Fig. IS6b . The two dashed lines in Fig. IS6b correspond to the factors XlC^i+^p) (left) and X2(^2 — W P ) 
(right). They have small values at off-resonance but cannot be neglected for stronger pump amplitude because the 
transition amplitude is proportional to V p 2 . The detailed calculation will be reported elsewhere. The series is given 
by 

xx(u) = irFU n 8{u - u d ) ^}^ 2) + ("d -► -w d ), (S28) 

where 

G ,(2) = Xi( w d)x2(wd + 2uj p ) [Axi(vd + u p )Ti + T 2 X2(^d + w p )A] 2 /16. 
We can similarly obtained the mode splitting for the second-order process as 

|Axi(w«i + w p )ri + r 2 X2(^d + w p )A| 



2Suj [ 



A| - Ti/wi + r 2 /w2| 

4((jJ2 — Wi)- s /WiCJ2 



(S29) 



The mode splitting for the third- and fourth-order processes can also be obtained as in the same way to the first- and 
second-order processes. Here, we only provide their final forms (the detailed calculation will be reported elsewhere): 

2fe( 3) „ 9Aj2Tj/qj + 2r 2 /^ 2 - (41^2 + A 2 )/^ 2 | 

128(w 2 — Wl) 2 - v /EJlCJ2 



(4) A| - 3rj/^3 + 9r 2 r 2 A^ 2 + 4A 2 iy^ 2 grirj/a^| - 4A 2 r 2 M^ 2 2 + 3rf/a|| 

36(w 2 - ^l) v^i^2 

If we assume that the piezoelectric detuning coefficient as a/2f2i = 2tt x 0.69 kHz/V, which is 8% smaller than that 
obtained when DC voltage is applied (Fig. IS2a). then we obtain the following theoretical values, which show very 
good agreement with the experiments as shown in Fig. 4 in the main text, 

2SJ 1] ~ 90.06 Vp_ p , 2SJ 2) ~ 60.23 V p 2 _ p , 2SJ 3) ~ 43.19 V^ p , 28u {i) ~ 31.43 V p _ p . 

The 8% difference in the detuning coefficient might be due to the aged shift or the voltage reduction due to the 
inductance of the electrical connection. 



Simulations of coherent mechanical oscillations 



The coherent mechanical oscillations can be simulated by directly solving time-dependent Eqs. IS12al and lS12bl The 
simulation results for the first- and second-order coupling are shown in Fig. IS7I with the experimental results. The 
experiments show very good agreement with the simulations. 




FIG. S7: Comparison of the pump voltage dependence of the coherent mechanical oscillations between the simulation and the 
experiment for the first- and second-order coupling. 



